A MICRO-MACRO PARAREAL ALGORITHM: APPLICATION TO 
SINGULARLY PERTURBED ORDINARY DIFFERENTIAL 

EQUATIONS 
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Abstract. We introduce a micro-macro parareal algorithm for the time-parallel integration of 
multiscale-in-time systems. The algorithm first computes a cheap, but inaccurate, solution using a 
coarse propagator (simulating an approximate slow macroscopic model), which is iteratively corrected 
using a fine-scale propagator (accurately simulating the full microscopic dynamics). This correction 
is done in parallel over many subintervals, thereby reducing the wall-clock time needed to obtain the 
solution, compared to the integration of the full microscopic model over the complete time interval. 
We provide a numerical analysis of the algorithm for a prototypical example of a micro-macro model, 
namely singularly perturbed ordinary differential equations. We show that the computed solution 
are better and better approximations of the full microscopic solution (when the parareal iterations 
proceed) only if special care is taken during the coupling of the microscopic and macroscopic levels of 
description. The error bound depends on the modeling error of the approximate macroscopic model. 
We illustrate these results with numerical experiments. 

1. Introduction. In many applications, a system is modeled using a high- 
dimensional system of differential equations that captures phenomena occurring at 
multiple time scales. Unfortunately, the computational cost of simulating such fine- 
scale systems (which we call microscopic in this work) on macroscopic time intervals 
is prohibitive, and one often resorts to low-dimensional, coarse-grained, effective mod- 
els (which we call macroscopic), in which the fast degrees of freedom are eliminated. 
Many methods have been proposed to obtain such macroscopic models, either ana- 
lytically (see e.g. [31] for a recent overview) or numerically. We refer, for instance, 
to the work on equation- free [3H or heterogeneous multiscale methods [5J [5] , and 
references therein. However, by construction, these macroscopic models only capture 
the original full microscopic dynamics approximately. 

Here, we present and analyze a numerical multiscale method that aims at ef- 
ficiently simulating the full microscopic dynamics (and not a macroscopic approxi- 
mation of it) over long time intervals, using an effective (approximate) macroscopic 
model as a predictor and the microscopic model as a corrector. To this end, we pro- 
pose a micro-macro version of the parareal algorithm [26) . The parareal algorithm was 
originally proposed to solve time-dependent problems using computations in parallel, 
aiming at exploiting the presence of multiple processors to reduce the real (wall-clock) 
time needed to obtain a solution on a long time interval. It is based on a decompo- 
sition of the time interval into subintervals, and makes use of a predictor-corrector 
strategy, in which the calculation of the corrections is performed concurrently on the 
different processors that are available. In what follows, we propose a version of this 
algorithm well-adapted to our multiscale-in-time context. 

For the sake of clarity, and to better describe our aim, we now present the parareal 
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algorithm in some detail. To fix the ideas, assume that the problem at hand is 



^ = /(«), «(0) = uo, te[0,T], (l.i) 

the exact flow of which is denoted u(t) = £t{uo). Suppose that we have at hand two 



propagators to integrate (1.1 ), J- At and GAt- The propagator J- At is a fine, expensive 
propagator, which accurately approximates the exact flow £At over the time range 
At, whereas the propagator GAt is a coarse propagator, which is a less accurate 
approximation of the exact flow. In turn, GAt is less expensive to simulate than T At- 



For example, J 7 At and Gai may correspond to integrating (1.1) over the time range 
At with a given discretization scheme, using either a small time step (for J" At) or a 
large time step (for GAt)- The parareal algorithm iteratively constructs a sequence 
of iV-tuplcs Uk = { u fe}i<„<jv ( w hh N = T/At), such that, at every iteration k > 0, 
it£ is an approximation of u(nAt). For k = 0, the initial approximation is obtained 
using the coarse propagator GAt- 

<±o = £At«L ), "°=o = wo- 
In the subsequent parareal iterations, the approximation is corrected using 



u 



n+l 
fc+1 



QmH +x ) + JAtK) - GAt{<), (1-2) 



with the initial condition u® k+1 = uq. The solution to ( |1.2[ ) can be very efficiently 
computed using the following procedure. Once the solution at parareal iteration k 
has been computed, we first compute the corrections FAt(uZ) ~ GAt(u^) in parallel 
over each subinterval [nAt, (n + l)At], < n < N — 1. We then only need to 
propagate these corrections sequentially, by adding GAt(u k+1 ) to the stored correction 
FAt(uk) ~ ^At(ufc)- This yields the solution at parareal iteration fe + 1. 

It has been shown (see e.g. [TJ [H [3J [5^) that, when k goes to infinity, the 
parareal solution converges to the reference solution, namely the solution given by 
the fine-scale propagator J- At used in a sequential fashion from the initial condition: 

Vn, < n < T/At, lim u£ = F$(uo). (1.3) 

k— foo 



The computational gain of the parareal algorithm stems from the fact that, in (1.2), 
the accurate simulations (using the fine-scale propagator IF At) are decoupled one from 
each other, and can therefore be executed in parallel on different processors. Suppose 
that the cost of a single evaluation of FAt is much larger than the cost of propagating 
the system according to GAt over the complete time range [0,T]. Assuming the cost 
of the fine-scale propagator FAt to be proportional to At, the cost of K iterations of 
the parareal algorithm is proportional to KAt. This cost is to be compared to the 
cost of computing the reference solution using the fine-scale propagator sequentially, 
which is proportional to NAt. The computational speed-up is thus N/K, which is 



larger than one if the number K of parareal iterations to obtain convergence in (1.3| 
is small enough. 

In this article, we propose and analyze a micro-macro version of the parareal 
algorithm. We assume that the variables in the microscopic model can be split into 
slow and fast components, and that we have at hand an approximate macroscopic 
model for the slow components under some time scale separation assumption (see 
Section [2] for the precise model we consider here) . In this setting, we will use the 
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parareal algorithm where the fine-scale propagator T&t is an integrator for the high- 
dimensional microscopic model, whereas the coarse propagator, here denoted Ca*, is 
an integrator of the low- dimensional, approximate macroscopic model (we use the 
notation Ca* rather than Q^ t to emphasize the fact that our coarse integrator acts 
on a system of smaller dimension than the reference one). The novelty therefore is 
to simultaneously use two models at different levels of description, rather than two 
discretizations of the same model. The cost of the coarse propagator is typically neg- 
ligible for two reasons: (i) the macroscopic model only contains the slow components 
of the evolution, and therefore allows for a larger time step; and (ii) the macroscopic 
model is low-dimensional, and therefore requires less work per time step. Again, the 
aim of the micro-macro parareal method is to speed up the computations (compared 
to a full microscopic simulation) by allowing the microscopic simulations starting from 
the different intermediate time instances nAt to be performed in parallel over each 
subinterval [nAt, (n + l)Ai], with < n < N - 1. 

As a model problem, we take the setting of singularly perturbed systems of ODEs. 
Such a model problem is a widely accepted first test case when proposing algorithms 
for problems with time-scale separation, see e.g. [5D]. We perform a numerical analysis 
of the algorithm we propose in a linear setting (see Section [2] for the description of the 
model problem, and Section [4] for the numerical analysis), and illustrate these results 
by numerical simulations in Section [5] However, our algorithm is not restricted to the 
linear setting, and we numerically observe in Section [6] that it performs equally well 
on a nonlinear test-case. 

Since its introduction in [26] . the parareal strategy has been applied to a wide 
range of problems, including fluid-structure interaction Navier-Stokes equation 
simulation |12j . reservoir simulation [15] , etc. The algorithm has been further analyzed 
in [3H 131] • Its stability has been investigated in [3 , 35J . An alternative formulation of 
the algorithm has been proposed in 0], or, equivalently, in pQ in a simplified setting. 
We refer to [H] for a reformulation in a more general setting that relates the parareal 
strategy to earlier time-parallel algorithms, such as multiple shooting (see e.g. [2"Tll33) ) 
or multigrid waveform relaxation (see e.g. [27j 36 ) approaches. Several variants of the 
algorithm have been proposed, for instance in [IT] [16] (see also [2] in the context of 
stochastic differential equations). The numerical analysis of the algorithm has been 
first performed for linear initial-value problems. A numerical analysis in a nonlinear 
context has been proposed in [13) . 

A micro-macro version of the parareal algorithm, similar to what is presented in 
this article, has already been considered in a number of works. The authors of [5j 128) 
consider a singularly perturbed system of ordinary differential equations (ODEs) at 
the microscopic level and the limiting differential-algebraic equation (DAE) at the 
macroscopic level. In these two works, the coarse propagator contains all degrees 
of freedom in the system. The slow degrees of freedom are evolved according to 
a differential equation, and the fast degrees of freedom are evolved using algebraic 
constraints (they somehow instantaneously adapt to the values of the slow degrees of 
freedom). In contrast, our approach completely eliminates the fast variables from the 
coarse propagator, and only evolves the slow variables. This difference has a number 
of consequences: 

• The coarse propagator in the algorithms proposed here is cheaper than that 
of jS] [25] (because it contains less degrees of freedom) and more convenient 
(because the coarse propagator simulates an ODE rather than a DAE); 

• The algorithms proposed here require operators to reconstruct microscopic 
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states from macroscopic ones, while the algorithm in (5] 128] can simply use 



the parareal iteration (1.2). This also influences the convergence behavior. 
A detailed comparison between our algorithms and that proposed in 5 ( 28 is given 
in Section I5~3l 

Other micro-macro parareal algorithms have also been proposed, in contexts dif- 
ferent from ours. In 10J , a parareal algorithm for multiscale stochastic chemical 
kinetics is presented, in which the macroscopic level uses the mean-field limiting 
ODE. In [32], the parareal algorithm is used with a kinetic Monte Carlo model at 
the macroscopic level and molecular dynamics at the microscopic level. 

Our article is organized as follows. In Section [2j we present the singularly per- 
turbed ODE that is considered here as a model problem, and state some bounds on its 
solution (The proof of these bounds is postponed until Appendix |A]) . Subsequently, 
in Section [3j we introduce two micro-macro parareal algorithms. The coupling be- 
tween the microscopic and macroscopic levels of description is done using a restriction 
operator (to go from the microscopic to the macroscopic level), and either a lifting 
(Algorithm 1) or a matching (Algorithm 2) operator (to go from the macroscopic 
to the microscopic level). This coupling ensures that the numerical solution remains 
consistent across both levels of description (see Section [3]) . The two algorithms we 
introduce in Section [3. 2| only differ in how the levels of description are coupled to each 
other. Algorithm [I] will turn out to be inaccurate, whereas Algorithm [2] is extremely 



accurate. For the sake of comparison, we discuss in Section 3.3 the scheme proposed 
in [3 [28], that we denote here Algorithm 3. Section [4] contains a detailed numerical 
analysis of these three algorithms, when applied to the linear model problem presented 
in Section [2] and when the dynamics at both microscopic and macroscopic levels of 
description are exactly integrated. This setting enlightens the effect of how the two 
levels of description are coupled on the convergence of the algorithms. We show how 
the modeling error of the approximate macroscopic model affects the accuracy. In 
particular, the micro-macro parareal algorithm we introduce is a precise approxima- 
tion of the full microscopic solution only if special care is taken during the coupling 
of the microscopic and macroscopic levels of description, as is done in Algorithm [2] 
The analysis is illustrated by numerical experiments in Section [5] where, in addition, 
we numerically investigate the effect of time discretization. Some numerical results 
on nonlinear problems are presented in Section [6] We observe there the same good 
properties of Algorithm [2] as on linear problems. We conclude in Section [7] with some 
final remarks and a discussion of possible future research. 

2. Model problem. In this section, we describe the microscopic model problem 
considered in this work, as well as its macroscopic limit. 
Consider the dynamics 

x = ax + p T y, y=-(qx- Ay) , (2.1) 
e 

where i£l and y G M d_1 are the state variables, and a G K, p G R d_1 , q G R^ 1 and 
A G M^ -1 )*^ -1 ) are parameters. This dynamics models the evolution of a system 
described by the state variable u — (x, y) G M. d , where the slow and fast components 
are x and y, respectively. We denote the initial condition by u(0) = (x(0),y(0)) = 
(xo,yo) = uq. The dynamics can be compactly written as 



u = B^u, 
4 



(2.2) 



where 



B e = 



a 



V 1 

-A/e 



In the following, we assume that the fast component of the system has a simple 
dissipative structure: 



We assume A to be a matrix with eigenvalues A,; £ C (1 < i < d - 1) 
satisfying Re(Ai) > A_ for any 1 < i < d — 1, for some A_ > 0. 



(2.3) 



Under this assumption, for each fixed value x = x* of the slow component, the 
dynamics of y, obeying the equation 

1 



satisfies 



y = - {qx* - Ay) 



lim y(t) = (A~ 1 q) x* . 



The dynamics of the fast component y, for fixed slow component x = x* , is thus 
exponentially stable for all x* . It is then known (see Lemma[2]below and, for example, 
[31] and references therein) that, in the limit e goes to zero, the solution x(t) to (2.1 1 
converges, on finite time intervals, to the solution X(t) of 



X = XX, X(0) = x , A 



■p A q. 



(2.4) 



Comparing (2.1 1 with ( 2.4 ) , one can see that the microscopic time-scale (namely the 
typical time-step required to integrate the full microscopic dynamics ( |2.1[ )) is of the 
order of e, whereas the macroscopic time-scale (namely the typical time-step required 
to integrate the approximate macroscopic dynamics (2.4|) is independent of e. 

Remark 1. The asymptotic result that we mentioned above on the system (2.1 1 
holds for more general cases. For instance, consider the dynamics 

1 



x = f{x, y), y = - (r)(x) - Ay) 



with again x G M, y G 



pd-l 



and A E R^- 1 )*^- 1 ) , and where f : M x 



pd-l 



(2.5) 
I and 



td-l 



are two given, possibly nonlinear functions. Under Assumption (2.3) , 



T] : K 

the solution x(t) to (2.5) converges to X(t), solution to 

X = F(X), X(0) = x , F(X) = f(X, A-^iX)). 

This result can also be extended to more general nonlinear cases \34V - 
For future reference, we introduce the exact time evolution operators, 

u(t* + At) = $ At (u(t*)) , X{t* + At) = PAt (X(t*)) , 



corresponding to (2.2) and (2.4), respectively. These equations are linear, hence the 
operators ^a* and pAt are linear: 



$At = exp(B £ Ai) e 
PAt = exp(XAt) e K 



[dxd 



(2.6) 
(2.7) 
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We now state some bounds on the solutions of (2.1), that will be useful in Sec- 
tion]^ when proving error bounds on the algorithms we propose. 

Lemma 2. Consider the linear system (2.1) over the time range [0, T], with the 
initial condition x(0) = Xq, y(0) = yo- Introduce z(t) — y(t) — A~ l qx(t) £ 1 and 
zq = z(0). Under Assumption ( |2.3| , there exist eq € (0, 1) and C > 0, that both only 
depend on A, q, p, a and T , such that, for all e < eq, 

sup \x{t) - x Q exp(A*)| < Ce(|sB | + \\zq\\), (2.8) 
te[o,T] 

sup \\z(t)-eiq>(-At/e)z \\<Ce(\x Q \ + \\zo\\). (2.9) 

te[o,T] 



Set 



Then, for all e < 6q, we have 



tf L = ^Hl/e). (2.10) 



sup \\z(t)\\<Ce(\x \ + \\z \\). (2.11) 

Hence, up to a boundary layer of size tf h , z(t) is of order e, and the state u(t) of the 
system is at a distance of the order of e of the manifold 

S := {u = (x,y) eR d ; y = A^qx] . (2.12) 



We call the manifold E the slow manifold. Note that the bound (2.11) is sharp 



in the sense that, after the initial time boundary layer, z(t) is of order e and not 
smaller. This can be checked for example on the analytically solvable system x — —x, 

y = (x- y)/e. 

An important consequence of the above lemma is that the microscopic solution 
u{t) = (x(t), y(t)) remains bounded, independently of e, on the time range [0, T\. The 
following result, which will be used repeatedly in the sequel, follows immediately from 
Lemma [2} 

Corollary 3. Consider the linear system (2.1) over the time range [0, T], with 
initial condition ir(0) = Xq, y(0) = yo. Under Assumption (2.3), there exist £q € (0, 1) 



and C > 0, that both only depend on A, q, p, a and T, such that, for all e < eo, we 
have 



sup \x{t)\ < C (\x \ + e\\y a \\) , 

t£[0,T] 

sup \\y{t)\\ < C (\x \ + e\\y \\) , 
te[tf L ,T] 



(2.13) 
(2.14) 



where the size tf of the boundary layer is defined by (2.10) 



The proofs of these standard results are postponed until Appendix [A] In view 



of (2.8), we see that, in the limit when e goes to zero, the macroscopic dynamics (2.4) 



is exact. The aim of the algorithms we investigate below is to use these macroscopic 



dynamics to speed up the computation of the solution of the original model (2.1 ), for 
a fixed small but non-zero value of e. 
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3. Micro-macro parareal algorithms. In this section, we describe two micro- 
macro parareal algorithms. As will become clear from the analysis in the forthcoming 
sections, the first one based on a lifting operator is inaccurate, whereas the second 
one based on a matching operator is extremely accurate. Both are generalizations 
of the parareal algorithm proposed in |26j . Our formulation follows most closely the 
description in pQ. We first introduce the necessary notation in Section 3.1 and we 



subsequently outline both algorithms in Section |3.2| For the sake of comparison, we 



also discuss in Section 3.3 the scheme proposed in [5JI2E]- Let us emphasize that the 
two algorithms we introduce are not restricted to the linear system (2.1), and apply 
to any system of the form 



x = f(x,y), y = -g(x,y), 

where x £ M. s is a slow component (s £ N*), y £ R m is a fast component (m £ N*), 
and where the associated macroscopic dynamics (obtained in the limit of infinite time 
scale separation between the slow and the fast components, namely in the limit when 
e goes to zero) reads X = F(X). 

3.1. Notation. We introduce a time discretization (t n )^ =0 , with t n = nAt and 
NAt — T. Let u n — (x n , y n ) ps u{t n ) — (x(t n ), y(t n )) be the numerical approximation 
of the solution of the microscopic model (|2.1[), and let X n sa X(t n ) be that of the 



solution to the macroscopic model (2.4) 



Fine-scale and coarse propagators. The micro-macro parareal algorithm makes 
use of two propagators. First, we need a fine-scale propagator, that advances the 



microscopic model (2.1) over a time-range At 



•FAt(« n )- 



(3.1) 



To perform this, we may consider that we have at hand the exact propagator of the 



equation (2.1 ), in which case JF^t = $At, where $At is defined by (2.6). Alternatively, 



we may resort to a numerical integration of the dynamics (2.1) (using for example 
forward or backward Euler discretizations) over the time range At, using several steps 
of size St. Typically, in the context of a system like (2.1 ), one would need St to be of 
the order of e to obtain accurate results. 



Second, we need a coarse propagator for the macroscopic model (2.4) 

X n+1 =C At (X n ), 



(3.2) 



where again we may assume that we can exactly integrate (2.4) and hence choose 
CAt = pAtj see (2.7). Alternatively, one may resort to a numerical integration of the 



dynamics (2.4), for which we can use a time-step independent of e to obtain accurate 
results. 

Restriction, lifting and matching operators. The parareal algorithm iteratively 
uses the fine-scale and the coarse propagators. In this work, these two propagators 
correspond to different descriptions of the system, either microscopic (using u £ K d ) 
or macroscopic (using X £ K.) . We thus need a way to go from one description to the 
other, as we discuss now. 

We first introduce the restriction operator 



K 



R d 



I — ^ X, 



which maps a microscopic state to the corresponding macroscopic state. For nota- 
tional convenience, we also introduce the complement of the restriction operator, 

n± ( R d -> R^ 1 

' \ u= (x,y) H> y, 

such that we can write u = (x, y) = (TZu, TZ^u). 

Conversely, we will also need to reconstruct a microscopic state from a given 
macroscopic state. In contrast to the restriction operator, there is no unique obvious 
way to define this operator. We introduce two such operators, a lifting operator and 
a matching operator. 

Definition 4. A lifting operator C is an operator 

( R -)• R d 
■\ X h> u = C{X) 

that creates a microscopic state that is uniquely determined by a given macroscopic 
state and satisfies the consistency property 

KoC = l&. (3.3) 



A possible choice is to take C(X) such that 

TZ(C{X)) = X and C{X) G E, (3.4) 

where E is the slow manifold associated to the multiscale problem. 

In connection with the system (2.1 ), an example (and this is the choice we make 
in this work) is to choose 



C(X) = (X, (A- 1 q)X). 



(3.5) 



This choice indeed satisfies (3.3) and ( |3.4[ ), where the slow manifold E of the sys- 
tem |2T]) is defined by ( |2.12[ ). 

Remark 5. Other lifting operators can be introduced, using for example the 
constrained runs algorithm [18]. As soon as the lifting operator C is specified, u is 
uniquely determined by X: the lifting operator enforces a closure approximation on 
the microscopic state. 

Alternatively, one may reconstruct a microscopic state using a matching operator. 

Definition 6. A matching operator is an operator 



V 



(X,v) ^ V x (v), 



that satisfies 



X = (K o V) {X, v) for any v e R d and X G 



(3.6) 



and 



Vu G R d such that TZ(u) = X, V x (u) = u, 
8 



or, equivalently, 



V(K(u),u) = u. 



(3.7) 



In contrast with a lifting operator, a matching operator requires a microscopic state 
as an input, and not only a macroscopic state. 



The consistency property (3.6) may be seen as the equivalent for V of the prop- 



erty (3.3) for C. We also note that, in view of (3.7), a microscopic state u which 



is already consistent with the macroscopic value X is unaltered by the operator 
Vx- Combining (3.6) and (3.7), we observe that Vx ° Vx = Vx'- the operator 
Vx '■ K d —> K d is thus a projection operator onto microscopic states u € R d that 
satisfy lZ(u) — X. One may thus think of Vx as a projection operator that projects 
a microscopic state v to a microscopic state u = Vx(v), such that lZ(u) — X and u 
is as close to v as possible, in a sense to be made precise for the problem at hand. 

In the following, we require in addition the following continuity property on V: 
there exists C > such that, for all X g R, Y e R, u € R d and v G K d , 



\\V(X,u)-V(Y,v)\\<C \\u-v\\ + \X -Y 



(3.8) 



For the analysis of the algorithms described below, we only require V to sat- 
isfy (3.6), (3.7) and (3.8), and do not make any additional assumptions (see Section|4]). 
For the numerical experiments reported on in Section [5j we choose, in the context of 
the system (2.1 ), 



V x (v) := (X,K x v), 



(3.9) 



which consists in keeping the fast variables from v, while imposing the slow variable 
to be equal to X. This choice fulfills all the above conditions (3.6), (3.7) and (3.8). 



Remark 7. The term matching operator has been chosen in reminiscence of the 
term "moment matching" that is commonly used in the Monte Carlo community, see 
e.g. J®. 

3.2. Algorithms 1 and 2. The parareal algorithm iteratively constructs ap- 
proximations on the whole time domain [0, T]. We denote by u%, resp. XJ} 7 the 
approximate microscopic, resp. macroscopic, solution at time t n , obtained at the fc-th 
parareal iteration. 

The first algorithm we consider is the following. 

Algorithm 1. Let u(0) = u be the initial condition. 
1. Initialization: 

a) Compute {Xfi } < n <jv sequentially by using the coarse propagator: 



X$=1l(u ), 



A™ +1 = C 



At 



TO- 



b) Lift the macroscopic approximation to the microscopic level: 



li 



and, for alll<n<N, u% = C(X$). 



are 



2. Assume that, for some k > 0, the sequences {u^} 0<n<N and {X£ } 0<n<N 
known. Compute these sequences at the iteration k + 1 by the following steps: 
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a) For all < n < N — 1, compute (in parallel) using the coarse and the 
fine-scale propagators 



Xl +1 =C M (X2), ^ +1 =7- At «). (3.10) 

b ) For all < n < N — X, evaluate the jumps ( the difference between the two 
propagated values) at the macroscopic level: 

j^ = n{ui +1 )-x n k +1 . (3.ii) 

c) Compute {X k+ i} Q<n<N sequentially by 

X° k+1 = K(u ), X'l+l = C At (X% +1 ) + J» +1 . (3.12) 

d) Compute { u k+i} o<n<N 1 °V lifting the macroscopic solution: 

u a k+1 =u and, for all < n < N - I, u' k l +l = C(X k l +l). (3.13) 
We can recast the above algorithm as 

KXl = c{c At (11 + n (T At «)) - C At (K «)) ) . (3.14) 



Notice that this cannot be recast in the form of the original parareal algorithm (1.2). 
The above algorithm uses the following paradigm: each time we need to reconstruct 
a full microscopic solution u from a given macroscopic state X, we use the lifting 



operator C. For example, for the system (2.1) and £ given by J3.5|, this amounts to 



by |3j} , 
2.12|. 



creating a microscopic state exactly on the slow manifold ( 2.12[ ). 

We will see in the sequel that this algorithm leads to disappointing results. In 
particular, Algorithm [T] does not retain one of the properties of the parareal algorithm 
as originally proposed in |26j . namely that the numerical trajectory is exact on the 
first k subintcrvals in time after k iterations of the parareal algorithm. 

A much better algorithm is the following: 

Algorithm 2. Let u(0) = u be the initial condition. 

1. Initialization: proceed as in Step 1 of Algorithm^ 

2. Assume that, for some k > 0, the sequences { u k}o< n <N an< ^ ^-X k } Q<n<N are 
known. To compute these sequences at the iteration k + 1, 

• Proceed as in Steps 2a, 2b and 2c of Algorithm^ 

• Compute {u k Xl} 0<n<N x by matching the result of the local microscopic 
computation, u k +1 , on the corrected macroscopic state X k +1 ■ 



L fc+i • 



T n+1\ 



ui +1 = u and, for all < n < N - 1 , = V(X r k 1 ^ , u k 

(3-15) 

The only difference between Algorithms [T] and [2] is how we reconstruct the micro- 
scopic solution u k Xl ■ In Algorithm [IJ we simply choose wjj+i on the slow manifold 
defined by X^+l (see pl3| >). In Algorithm B we use the quantity u k +1 , which is the 



end point of a microscopic trajectory between times nAt and (n + l)Ai, and match 
this state onto the corrected macroscopic state X^\: obtained at the latest parareal 
iteration. 
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At the initial iteration k = 0, since no microscopic computation has been done, 
we cannot use the matching operator V to reconstruct a fine-scale solution. We thus 
resort to the lifting operator C. 

Algorithm [2] can be recast as 

ultl = v(c At (TZ + Tl (J- At «)) - C At (TZ «)) , 7-AtK)) , (3.16) 



which is to be compared with the original parareal algorithm (1.2 1 and (3.14) for 
Algorithm [T] For the linear system (2.1) and the choice (3.9) of matching operator, 
the equation ( |3. 16 ) can be further simplified to 



i£l = -F At «) + (1, 0) T (c At (TZ (ut +1 )) C At (TZ «)) ) . (3.17) 



This is exactly (1.2) with T A t as the fine propagator and (1, 0) T C A tTZ as the coarse 
propagator. 

Note that, in view of (3.3) and ( |3.13 ) (respectively (3.6) and (3.15)), the trajec- 
tories computed using Algorithm [T] (respectively Algorithm |2| satisfy 



Vfc > 0, Vn > 0, XL 



1l(ut). 



(3.18) 



At any parareal iteration fc, the macroscopic trajectory is consistent with the micro- 
scopic trajectory. 

3.3. Comparison of Algorithms[l]and[2]with that of [5ll28] . As underlined 
in the introduction, a micro-macro version of the parareal algorithm has already been 
proposed in [28j [5] . In these works, the coarse propagator is an integrator of a reduced 
(DAE) model that contains all degrees of freedom in the system (both the fast and 
slow ones), in contrast to our algorithms, where the coarse propagator is an integrator 
for the effective dynamics of the slow degrees of freedom. 

For the model problem (2.1), the reduced DAE considered in [251 \E\ takes the 
form 



T 

P V: 



Ay = qx. 



(3.19) 



The coarse propagator of [28j [5] is an integrator Q At of (3.19). This coarse inte 



grator is combined with a fine-scale integrator J- At of (2.1) in the parareal fashion, 
following (TOl). 



The obtained scheme, that we denote here Algorithm 3, differs from our Algo- 
rithm [2] in its treatment of the fast degrees of freedom. To show this, we note that, 
specifically for the model problem (2.1), an exact propagator for (3.19) can be ob- 



tained by first solving (2.4) exactly, and second solving the algebraic equation for y. 
Hence, we have 



G A t (u) — Co p At o TZ u = C o C A t olZu 



(3.20) 



where, we recall, C At is the coarse propagator used in Algorithms [T] and [2] Using ( 1.2 ) 
we write Algorithm 3 as follows: 



CI = ^tK) + g At {u n k+1 ) - g At «) 

= F At (ul) + C(C AI {TZ - C At (TZ «)) ) , 



(3.21) 
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which can be compared with (3.171 and with (3.14). Notice in particular that Algo- 
rithm [2] differs from Algorithm 3 in the choice of the coarse propagator. 



The three algorithms only differ in how the microscopic and macroscopic levels 
of description are coupled in the parareal iterations. These differences, however, have 
implications on (i) the computational complexity of the methods; (ii) the way they 
generalize to more complex multiscale systems; and (iii) the convergence behavior. 
The convergence properties of the three algorithms are analyzed in Section [4j We here 
briefly comment on the other two aspects. First, the computational complexity of the 
coarse propagator in Algorithm 3 is significantly higher than that of Algorithms [l] 
and [2] due to the presence of the fast degrees of freedom, which requires solving a 
large linear system in addition to the time-stepping of the slow degrees of freedom. 

Second, in more complex situations, for instance when the microscopic and macro- 
scopic systems arc nonlinear, Algorithm 3 may require the use of a time integrator 
for DAEs. Although many such integrators exist, they are usually implicit, and less 
convenient than ODE solvers. In those cases, Algorithms [T] and [2] only require a 
reasonable model to propagate the macroscopic variables. In both cases, one may 
resort to computational multiscale methods that approximate the evolution of the 
approximate macroscopic model by using short microscopic simulations. The coarse 
propagator required in Algorithms [l] and [2] can be replaced by a coarse projective 
integration approach [221 H3] • The coarse propagator for the DAE system required 
in Algorithm 3 can be replaced by a projective integration method [T7J [TS]. Remark 
that the computational cost of both methods is not identical: projective integration 
requires a computational cost of 0(log(l/e)), whereas the computational cost of coarse 
projective integration is independent of e. This shows again that Algorithms [l] and [2] 
are cheaper to implement than Algorithm 3. We will see in the next Section to what 
extent the higher computational cost of Algorithm 3 allows for a better accuracy. 

4. Analysis. In this section, we analyze the convergence of the two micro-macro 



parareal algorithms introduced above on the linear model problem (2.1 1. We also give 



a detailed analysis for Algorithm 3, introduced in jS][2H]. To keep the analysis simple, 
we focus on the error due to the fact that the models are different at the macroscopic 
and microscopic levels. We thus track the dependency of the error bounds on the 
parameter e, and consider, at both levels, the exact propagators (2.6) and (2.7). 



Thus, the fine-scale and coarse propagators in (3.1) and (3.2) are given by 



^At(w) = ®Atu, C A t(X) = pAtX, 



for a fixed At, which is chosen typically much larger than e (so that At is a macroscopic 
time-scale). We recall that the lifting operator C is defined by (3.5 1, and that we work 



with a matching operator V satisfying (3.6), (3.7) and (3.8) 
We first derive an error recursion formula in Section |4.1 



Using this formula, 

we derive a sharp error bound on the trajectories computed by Algorithm [I] where 
the microscopic state is reconstructed using the lifting operator £ (see Section |4~2| ). 
We next turn to Algorithm [2j where the microscopic state is reconstructed using 
a matching operator V. We first show that, at a given parareal iteration k, the 
computed trajectories (both at the macro and the micro scales) are exact up to the 
time kAt (see Section 4.3.1), reproducing thereby a property of the standard parareal 
algorithm. We subsequently derive a sharp error bound in terms of e, showing that, 
at iteration k, Algorithm [2] converges to the exact so lution of the full microscopic 

for precise statements). 



system with an error of the order of e fc / 2 ( see Section 
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4.3.2 



These two properties (exactness of the trajectories up to time kAt after k iterations, 
and improvement of the convergence rate to the exact solution as k increases) are not 
satisfied for Algorithm[l] We eventually consider Algorithm 3. Being based on (1.2) 



this algorithm automatically satisfies the local exactness property. We then prove a 
sharp error bound in terms of e, showing, in agreement with [28] . that, at iteration 
k, Algorithm 3 converges to the exa ct so lution of the full microscopic system with an 



error of the order of e k (see Section 4.4 for precise statements). 

The analysis below closely follows that of [26], but is significantly extended. We 
explicitly relate to the case considered in |26j when appropriate. 
Before proceeding, we introduce two notions of error: 

Definition 8 (Microscopic error). Let u{t n ) be the exact microscopic solution 



of (2.2) at time t„ = nAt, and let vll be the parareal microscopic solution after k 



parareal iterations, using Algorithm^ or[| The microscopic error 

et = ut- u(t n ) (4.1) 

is defined as the difference of the solutions at the microscopic level. 

Definition 9 (Macroscopic error). Let u(t n ) be the exact microscopic solution 



of (2.2) at time t n — nAt, and let X£ be the parareal macroscopic solution after k 



parareal iterations, using Algorithm^ or[l[ The macroscopic 

El = XI - Ku(t n ) (4.2) 
is defined as the difference of the solutions at the macroscopic level. 



Note that, in view of (3.18) and using the linearity of 1Z, we have 

E£=Ket. (4.3) 

4.1. Error recursion formula. A first step in the analysis of the algorithms 
described above is the derivation of a recursion formula for the error, which is valid for 



both algorithms and for any choice of the operators TZ, C and V . Starting from (3.12 ) 



and (3.11 1, we write X% +1 as a function of the microscopic and macroscopic solutions 



at the parareal iteration k: for n > 2, 

Xn n ( \rn— 1\ i jn 

k+1 — ^At(^ k +1 ) + Jk 



k 



= pAtX^l + (T^At^r 1 - PA* A 

= n^ At u^ + PAt (i^-ir 1 ) 

= ft* At ur x + pa* (n<s> At u n k - 2 + pa* (xi-l - xr 2 ) - xi- 1 ) 

1 , _ „, n ~2 xrn— 1\ i 2 / yn—2 \rn—2\ 



= ft^r 1 + pa* (^A^r 2 - xr 1 ) + P \ t (a^ 2 - n 

n-1 

= rc* A t«r 1 + e p p A t (^A^r*- 1 - a 



-n—p 



0=1 

-1 



■R^muT 1 + E pa7 (k*m«V - xi) . (4.4) 
P =i 



This formula is also valid for n = 1 using the convention ^ ■ = 0. Note that we 
have used the linearity of the coarse propagator. We then obtain a recursion for the 
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macroscopic error, using the linearity of the fine-scale propagator: 



n-l 



P =i 

n-l 

p=l 

n-l 



n— 2 n— 1 

p—0 p—1 

= 'j2p n At P - 1 mAte p k - PAt El), (4.5) 

p=l 

where, in the last line, we have used the fact that the microscopic error at t = 
vanishes: e° = for any fc. 



We remark that the formula (4.5) is not closed, in the sense that it couples 
the macroscopic error at parareal iteration k + 1 to both the macroscopic and the 
microscopic errors at parareal iteration k. To close the formula, and to transform it 
into specific bounds on the errors, we will need to make use of specific properties of 
<&At and of the lifting, matching and restriction operators C, V and 1Z. This is where 
the analysis of Algorithms [T] and [2] differ. 

Remark 10. Using ( |4.5[ ) ; it is possible to recover standard error bounds on the 
parareal algorithm, when the microscopic and the macroscopic models are linear and 
written at the same level of description, using a common state variable u G R, as 
in \26\/ for example. In this case, we have 1Z = C = V = Id (there is only one model, 
and one level of description), and ejj = E r k l . The coarse and fine- scale propagators are 
linear operators, denoted respectively byC At {u) = GAt{u) = P At u and J- At ( u ) — PAt u - 
Since u is scalar, the propagators are simply multiplications by two scalars p^ t and 
Pa*' The equation (4.5) then reads 



E% +1 = Y,(PAt) n - p - 1 (Pit - Pit) El (4.6) 

Assume, as in the classical analysis of the parareal algorithm presented in [26], that the 
fine-scale propagator is exact, whereas the coarse propagator is a scheme of order s: 
\p At — p%t\ — 0(At s+1 ). We consider a range of At such that p^ t > (which is 
possible since p^ t = 1 + O(At)). Fix a time range [0, T]. We show that, using ( |4.6| , 
one can recover the classical result of \2&j : at any parareal iteration k, there exists c k 
such that, for any At, 

sup \£%\ < c k At s( - k+1 l (4.7) 

0<nAt<T 

This bound is satisfied at k — since the coarse propagator is of order s. Assume now 



that (4.7) holds at some parareal iteration k. We then deduce from (4.6) that, for all 
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> such that nAt < T, 



n-l N-l 

\E^ +1 \ < c k At°( k +V \ P l t - p %\ ^(pltT-*- 1 < Cc k At^At^ J2(PAt) P - 

p—1 p—0 

(4.8) 

Remark now that 



N-l 

£(/>: 

p=0 



pit -i 



Since the fine-scale propagator is exact, we have (p& t ) N — PpjAt ~ Pt> which is 
independent of At. Thus 



N-l 
p=0 



CAt s + C C 

\PAt 



(4.9) 



Collecting (4.8) and (4.9), we deduce (4.7) at the parareal iteration k + 1. This 



concludes the proof. 

4.2. Error bounds for Algorithm [TJ We consider Algorithm [l] where the 
reconstruction at each parareal iteration is done using the lifting operator £ defined 
by (3.5). We show in this section that the accuracy of the numerical trajectory 



does not improve (neither at the microscale nor at the macroscale) as the number of 
parareal iterations k goes to infinity. 

Theorem 11. Consider Algorithm^ where J-a* is the exact propagator of the 
microscopic problem (2.1), Ca* is the exact propagator of the associated macroscopic 
problem (2.4), and C is the lifting operator defined by (3.5). We fix the time range 
[0, T], and recall that the size of the boundary layer tf h in (2.1) is defined by (2.10). 



Then, there exists £q £ (0, 1), that only depends on A, q, p, a and T, such that, 
for all e < €q and all At > tf h , there exists C, that depends on A, q, p, a, At and T 
such that 



sup \E£\<Ce and, for all k > 1, sup \E%\ < Ct\ 

0<n<N Q<n<N 



for all k > 0, sup \\v k 

0<n<N 



et\\ < Ce, 



(4.10) 
(4.11) 



where N = T/ At and where the macroscopic (resp. microscopic) error E^ (resp. e k ) 



is defined by (4.2) (resp. (4.1 )). Note that C is independent from e and k. 

The numerical experiments described in Section [5] show that these error estimates 
are sharp. Recall that if L = Celn(l/e) for some constant C only depending on the 
matrix A of (pO} (sec (|2.10[)). The assumption At > tf L = Celn(l/e) is therefore 



automatically satisfied for sufficiently small e, and in particular when the time-step 
At is of the order of the macroscopic time scale. 



Proof. Using the definitions (4.1) and (4.2) of the microscopic and the macro- 
scopic errors, and the fact that the microscopic state is reconstructed via the lifting 
operator C, 
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we have 

£ El = c xt ~ £K $l t u = u n k - £K $ n At u = e n k + $ n At u - £K $ n At u , 
or, equivalently, 

e n k =£E%- (Id - £11) <S> n At u . (4.12) 



As a consequence, we can write the recursion ( |4.5| ) for E^ +1 in terms of E v k only, by 
eliminating the microscopic errors e p k . We have 

TVf> At e p k - p At E p k = K$ At [£ E p ~ (Id -£TZ) § At u ] - p At E p 

= (R& At £ - p At ) E p - K<5> At (Id -OR) <S> p At u . (4.13) 



The first term in (4.131 stems from the difference between the macroscopic evolution 
of the microscopic system and the evolution of the approximate macroscopic equation. 
The second term stems from the difference in evolution between a microscopic state 
and the (unique) microscopic state that is obtained by lifting its restriction. To bound 
the first term in (4.13), we observe that 



\n$ At £-p At \ < Ce. 



Consider indeed the system (2.1) with the initial condition (xo,yo) = 
z a = (because £(xq) e S, see (3.5)), and we deduce from (|2.8[) that 



(4.14) 
C(xq). Then 



that reads 



\x(At) - x exp(XAt)\ < Ce\x \, 
\lZ$ At £(x ) ~ pAtXa\ < Ce\x \, 



from which we infer (4.14). 

We now turn to the second term of equation (4.13). We introduce the shorthand 
notation for the exact solution 

<f> p A t u = u p = (Ku p ,K x u p ) = {x p ,y p ) . 



First, using the definition (3.5) for £, we write 



(Id-£K)<S> p At u 







y p -(A- x q) x p 



Second, using ( 2.13[ ) with the initial condition xq = 0, y = y p — (A 1 q) x p , we get 

\K$ At {Id-£K)$ p At u \ <Ce \\y p - {A- X q)x p \\. (4.15) 

We are now left with bounding ||y p — (A~ 1 q)x p . We note that y p — (A~ 1 q)x p = 
z(pAt), thus, using (2.11) for the solution u(pAt) — & At uo, we deduce that 

\\y p - (A^q) x p \\ = \\z( P At)\\ < Ce (\x \ + \\z Q \\) < Ce. 

Note that we have used the fact that p > 1 and At > tf L , hence pAt > tf L . We then 
deduce from (4.15) that 



\K$ At (Id-£K)$ At u \ < Ce 2 
16 



(4.16) 



Collecting fl4.13| ), ( |4.14[ ) and ( |4.16[ ), we obtain 

|7e$ A ^-M^|<Ce(K|+e). 
where C only depends on A, q, p, a, T and At. 



Inserting this bound into the error recursion (4.5), and using that p/\ t > 



(see ( 2.7 )), we get 



\E^ +1 \<CeJ2pTr 1 (\E P k \+^) 



We now introduce E k '■= nrax <„<T/At \E k |, and write 

\E n k+1 \ < Ce (E k + e) £ p^ 1 = Ce + e) i 



PAt 



Let m 



0<n<N 1 — pAt 



which only depends on At, T and A. We obtain 



Ek+i < Cme 



(X + e) 



where Cm only depends on A, q, p, a, T and At (and is in particular independent of 
k and e). We thus have 

< E k < v k 

where the sequence {ffe} fcgN is recursively defined by v k+ i = Cme (v k + e) and vq = 
E , so that 

1 - (Cme) 4 



v k = E {) {Cme) k + Cme 2 



1 - Cme 



Note that the bound ( |2.8[ ) reads \x(t) - X(t)\ < Ce, hence v = E < Ce. 

Let us choose eo = l/(Cm). Notice that eo only depends on A, q, p, a, T and At. 
For any e £ (0, eo), the sequence v k has a limit as k goes to infinity and there exists 
C, independent of k and e, such that 

< En < Ce and Vfc > 1, < E k < v k < Ce 2 . 



This proves the bound (4.10) on the macroscopic error. 

To prove the error bound on the microscopic error, we notice, using the defini- 
tion (3.5 ) of £, that 

(Id-CTZ) $ At u = (Id-CTl) u{nAt) = (o,y(nAt) - (A~ 1 q)x(nAt)j . 

Since At > tf L , we deduce from (|2.11|) that 



Vn > 1, \\(ld-£K) $i t «o|| < Ce. 



(4.17) 



Collecting (4.121, (4.10) and (4.171, we obtain, for any k > 0, 



Vn > 1, 



<C\E%\ + \\(Id-Ol) $l t u \\ <Ce. 



Note that the microscopic error is always dominated by the lifting error (the second 
term of ( |4.12| ). 

Since, at any parareal iteration k, we start with the correct initial condition, we 
have el = and we thus have proved (4.11). □ 
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4.3. Error bounds for Algorithm [2} We now consider Algorithm [2j where 
the reconstruction at each parareal iteration is done using any matching operator V 
satisfying (3.6) and (3.7). The continuity assumption (3.8) will be added when needed. 
As pointed out above, we do not assume any specific expression for V here. We show 
in Section [4. 3. 2| that, in contrast to Algorithm [lj the convergence rate obtained with 
Algorithm [2] increases as the number of parareal iterations k increases. Before that, we 
show in Section |4 . 3 . 1 1 that , at a given parareal iteration fc, the computed trajectories 
(again both at the macro and the micro scales) are exact up to the time kAt. 

4.3.1. Local exactness of the algorithm. One of the important properties 
of the parareal algorithm (1.2) is that it results, after k parareal iterations, in a 
solution that is exact at all times up to kAt. The word "exact" here means that the 
parareal solution is equal to the solution that would have been obtained using only, 
in a sequential fashion, the fine-scale propagator up to time kAt. We now show that 
this exactness property holds for the micro-macro parareal algorithm we propose. 

Theorem 12. Consider Algorithm^ where T A t is the exact propagator of the 
microscopic problem (2.1), C A t is the exact propagator of the associated macroscopic 
problem (2.4), C is the lifting operator defined by (3.5) andV is a matching operator 
satisfying (3.6) and (3.7). 

Denote by vJJ, the microscopic solution obtained at the n-th time-step and k-th 
parareal iteration, using Algorithm^ Then, at any parareal iteration k > 1, we have 



\fp < k, 



$ At tt . 



(4.18) 



Proof. The proof goes by induction. Consider the parareal iteration k = 1. We 
obviously have = uq — $ A( uo- At time iteration n = 1, in view of (3.15), we have 



u\=V{Xl,ul), 



with (see ( |3.12[ )) 



■C At (X°)+K(ul) 



Cai(Xq) 



u\ = V{K{u\),ul) = u 1 = J"a*(uo) = $ AtU . 



n{ul). 



Hence, using the fundamental property (|3.7|), 



This proves (4.18) for k = 1. 

Assume now that, at some parareal iteration k > 1, we have (4.18). In view 
of (3.18), this implies that X^ = 7Z [$ At uo] for any p < k. Using (4.4) and the fact 
that $a*w? _ = ^At^o for a ll p < k, we deduce that 



Vn< + 
Hence, at the parareal iteration k 



T^A^r 1 



n<s>i t u . 



(4.19) 



1, the macroscopic solution is exact up to time 



(k + l)At. Using (3.15), we now write, for any n < k, 



where we have used ( 4.19[ ) and (3.10) in the first equality, the exactness assumption 
of the microscopic solution at iteration k in the second equality, and the fundamen- 
tal prop erty (|3.7[ ) of the matching operator V in the last equality. This proves the 
relation ( 4.18[ ) at the iteration k + 1 and concludes the proof. □ 

This result also directly follows from our above remark that, in its form (3.17), 
Algorithm 2 is of the form (1.2). 
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4.3.2. Error bounds. We now establish error bounds on Algorithm [2] that show 
that the microscopic solution converges towards the exact microscopic dynamics when 
the modeling error e decreases, and that the convergence rate improves as the number 
of parareal iterations k increases. This is in contrast with Algorithm [T] where the 
error does not improve even if k goes to infinity (see Section 4.2). With Algorithm K2| 



we recover the behavior of the standard parareal algorithm, as recalled in Remark n0\ 



(see e.g. (4.7)). 

Theorem 13. Consider Algorithm^ where J 7 ^ is the exact propagator of the 
microscopic problem (2.1), Ca* is the exact propagator of the associated macroscopic 



problem (2.4), C is the lifting operator defined by (3.5), andV is a matching operator 
satisfying (3.6), (3.7) and (|3.8[). We fix the time range [0,T], and recall that the size 



of the boundary layer tf in (2.1) is defined by (2.10) 



Then, there exists 6q € (0, 1), that only depends on A, q, p, a and T, such that, 
for all e < e and all At > ijp L , there exists a constant Cf., independent of e, such 
that 



for all k > 0, sup \E%\ < C k e 1+ ^ k/ ^ , 

0<n<N 

for all k > 0, sup ||eJJ|| < C fc e 1+Lfc/2J , 

0<n<N 



(4.20) 
(4.21) 



where N = T/ At and where the macroscopic (resp. microscopic) error Ej! (resp. e k ) 
is defined by (4.2) (resp. (4.1)). The constant Cfe is independent from e, but a priori 



depends on k, A, q, p, a, At and T . 



In (4.20) and (4.21), we used the notation: for any [~af| € Z and [^J G Z 

are respectively defined by: [^J < x < [x\ + 1 and \x~\ — 1 < x < \x] . 

The above result shows that the parareal iterations alternatingly improve the 
macroscopic and the microscopic errors by an order of magnitude in e. The numerical 
results of Section [5] show that (4.20) and (|4.21 ) are sharp error estimates. As already 
mentioned above, the assumption At > tf^ is automatically satisfied for sufficiently 
small e, in particular when the time-step At is of the order of the macroscopic time- 
scale. 



The bounds (4.20) and (4.21) show that, as k increases, the rate of convergence 
(with respect to e) of the error increases. The dependence of the constant C k in 
these two bounds on At and k will be analyzed in details on the numerical test case 



considered in Section 5.1 (see (5.4) and (5.5)) 



Proof. Using (3.15), (3.7) and the definition (4.1) of the microscopic error, we 



have 



Hence, using (3.8), we deduce that 



E 



fc+i 



< C $ 



Ate k 



E 



fe+i | 



Since At > tf L , we infer from ( |2.13| and ( |2.14[ ) that 



I^Aterl <C(\Ke 



-L„™-1 1 



C(\E, 



n— 1 1 
k 



-L„™-1 1 



(4.22) 
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where we have used (4.3). We then deduce from ( |4.22| that 
|K +1 || < C(\E^\ + e ||^Vr 1 || + \E% +1 \) < C(\E^\ +e\\e 



J k 



E 



k+l\) ■ 
(4.23) 

We now bound the macroscopic error E k+Il using the recursion formula (4.5), 
that reads 



(4.24) 



with T k : = lZ&Ate k — PAtE k . Consider the solution (x(t),y(t)) to the system (2.1) 
with initial condition u(0) = e p k , that is x(0) = K(e p k ) = E p k and y(0) = K^{e v k ). We 
then have 

T p = x(At)-X(At), 



where X(At) is the solution to (2.4) with initial condition X(0) = E p . In view 
of J2~8l, we have 



\T%\ = x(At) X(At) < Ce (\x(0)\ + 111/(0) - A^qx^W) < Ce (\E P \ + ||e£||) , 
where C is independent from p, k and e. We are now in position to use the recur- 



sion (4.24), from which we infer 



fell < E^a/" 1 1^1 ^ ^E^" 1 (Kl + KID- 

p— 1 J3— 1 



(4.25) 



We now prove the theorem by induction, using the two fundamental estimates ( 4.23[) 
and (4.25). At the parareal iteration k = 0, Algorithm [2] is identical to Algorithm [1] 
In view of (4.10) and (4.11), we thus have 

sup \Eq \ < C e and sup 1 1 eg" 1 1 < C e, 

0<n<N 0<n<N 



that is ( |4.20| and ( |4.21[ ) for k = 

Let us now assume that ( 4.20[ ) and (4.21) hold at any iteration k' < k, with k 
an even integer. We prove the bounds at iteration k + 1. Setting m = k/2 + 1 (so 
that \k/2\ = \k/2] = m- 1, |_(fc-l)/2j = m-2 and \(k- l)/2] = m- 1), we thus 
assume that 

sup \El_A <C fc _ l£ m , sup llejjl < Cfc-ie™- 1 , 

0<p<N 0<p<N 

sup |^| <C fc e m , sup K||<C fc e m . 

0<p<N Q<p<N 



Then, we infer from (4.25) that, for any < n < N, 

i _ n n— 1 ^t/^i ,m+l 
I pn I ^ ^i,^ ,m+l \ ^ n-p-1 ^ ^-.^ ,m+l 1 ^At s l^fc fc „ m+1 

p~l 1 - PAt 1 - PAt 

where Cfe+i is independent from e, but depends on k and At. We next deduce 
from ( |4.23[ ) that, for any < n < N, 

||e£ +1 || < C (C k e m + C k e m+1 + C k+1 e m+1 ) < C k+1 e m , 
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where C k +i is again independent from e, but depends on k. We thus have proved ( 4.20 1 
and (4.21 1 at iteration k + 1. 

We now assume that (4.20) and (4.21) hold at any iteration k' < k, with k an 
odd integer. We prove the bounds at iteration k + 1. Setting m = (k — l)/2 + 1 (so 
that lk/2\ =to-1, \k/2\ = m and [(k - l)/2j = |~(fc - l)/2] = m - 1), we thus 
assume that 



sup 

0<p<N 



sup \E p k \<C k e r ' 

0<p<N 



sup lie 

0<p<N 

sup 

0<p<N 



fc-1 



Using again equations (4.25) and (4.23), we find that, for any < n < N, 

m+1 and He^JI <C k+1 e m+ \ 



Pfc+1 



<c, 



'fe+1 1 



where C k +i is again independent from e. We thus have proved (4.20) and (4.21) at 
iteration k + 1. This concludes the proof. □ 

4.4. Error bounds for Algorithm 3. Since Algorithm 3 uses the standard 
parareal iteration (1.2), local exactness is automatically satisfied. We proceed to 
proving error bounds on Algorithm 3, which can be compared to those of Algorithm[2j 

Theorem 14. Consider Algorithm 3 given by (3.21), where J-'ai is the exact 
propagator of the microscopic problem (2.1), Q& t is the exact propagator of the asso- 
ciated macroscopic DAE problem (3.20), and C is the lifting operator defined by (3.5). 
We fix the time range [0, T], and recall that the size of the boundary layer tf L in (2.1 ) 
is defined by (2.10). 

Then, there exists cq £ (0, 1), that only depends on A, q, p, a and T , such that, 
for all e < £o and all At > if L , there exists a constant C k , independent of e, such 
that 



for all k>0, sup \E r k 

0<n<N 

for all k > 0, sup ||e£ 

0<n<N 



< C k e k+ \ 
<C k e k+ \ 



(4.26) 
(4.27) 



where N = T/ At and where the macroscopic (resp. microscopic) error E^ (resp. e r k l ) 
is defined by (4.2) (resp. (4.1)). The constant C k is independent from e, but a priori 
depends on k, A, q, p, a, At and T. 

These results are in agreement with [551 Theorem 2.1]. The numerical results of 
Section [5] show that (4.26) and (4.27) are sharp error estimates. 

The above result shows that, in contrast to Algorithm [2j Algorithm 3 improves 
the order of convergence (in terms of e) of both the macroscopic and the microscopic 
errors by an order of magnitude in e at each iteration. As noted in Section |3.3| 
this improved convergence rate comes at the price of a larger computational cost per 
iteration. 



Proof. Using (3.21 ) and the definition (4.1 ) of the microscopic error, we have 



-fc+i 



Using (4.3), we deduce from the above equation that 

-1 



-fc+i 



< 



I^K 



k+l 



(4.28) 
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The first term is decomposed as 



PAt l<-e k 



Since At > tf L , we infer an estimate on the first (resp. second) term of the above 
right-hand side using (2.8) (resp. ( 2.11[ )), resulting in 



I ^ n—X 



CpAtTle 



n-1 I 
k 



< Ce 



n-1 1 



(4.29) 



of e, such that 



Collecting (|4.28|) and (|4.29|), we deduce that there exists a constant C, independent 

<C(e\\el- l \\ + \E n k -{\). (4.30) 



-fc+i| 



This estimate is to be compared with (4.231 in the proof of Theorem 13 Using (4.30), 
the proof of Theorem 14 is completed via induction, in a way that is similar to (but 



simpler than) the proof of Theorem 13 □ 

5. Numerical experiments (linear test-case). In this section, we numeri- 
cally illustrate the above convergence results on a linear problem. We first consider 
the case when both the microscopic and the macroscopic models are integrated ex- 



actly (Section 5.1). We next consider the case when the macroscopic propagator is a 
forward Euler discretization, thus introducing some finite step-size error (Section 5.2). 
We consider the example system 




(5.1) 



which is of the form (2.1). The associated macroscopic, slow dynamics is given by ( 2.4 ) 
with A = — 1. The initial condition is x(0) = 1, yi(0) = 2/2(0) = 0, and we consider 
the solution on the interval [0, T] with T = NAt = 10. 

The fine-scale propagator J 7 At is the exact integrator of ( |5.1[ ). The coarse prop- 
agator Ca* is the exact integrator of (2.4) in Section 5.1 and a forward Euler dis- 



cretization of (2.4) in Section 5.2 We choose the parareal time-step At = 10 1 , and 
consider e e [10~ J , lO -1 ]. The lifting operator C is defined by (3.5), and we use the 



matching operator V defined by (3.9) 



We look at the relative macroscopic error \E£ 



scopic error ||e£/ || / ||u(T)|| at the final time T 
numbers k, satisfying < k < K. 



/ \x(T)\ and the relative micro- 
tjv = NAt for different iteration 



5.1. Results with exact microscopic and macroscopic integrations. In 

this section, we take both the fine-scale and the coarse propagators to be the exact 
integrators. 



5.1.1. Algorithm [T} We first consider Algorithm [T] (analyzed in Section 4.2), 
where the reconstruction at the end of each parareal iteration is done using the lifting 
operator C. 

We set the maximal number of parareal iterations at K = 2. Figure [5~T| shows 
the macroscopic and microscopic errors as a function of e for the chosen values of k. 
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We see that the macroscopic error is of the order of 0(e 2 ) as soon as k > 1 (and is of 
the order of 0(e) at k = 0). The macroscopic error at k — 2 is equal to that at k = 1. 
We also observe that the microscopic error is always of the order of e (for any fc), 
although the value of the error is smaller at k = 1 than at k = 0. These results are in 
agreement with Theorem 11 and confirm the fact that the accuracy of Algorithm [l] 
does not improve when k goes to infinity. 





Fig. 5.1. Algorithms/or the system \5.1\ , with exact fine-scale and coarse propagators: errors 
as a function of e for different values of k (left: macroscopic error; right: microscopic error). Note 
that the lines for k > 1 visually overlap. 



When e is too large, the macroscopic error is not anymore of the order of 0(e 2 ) 
at the iteration k > 1. This is due to the fact that the assumption At > tf h is no 
longer satisfied. Recall indeed that we keep At fixed, and tf h = Celn(l/e) increases 
if e increases. Hence, for too large values of e, the time step At is too small to correct 
for the initial boundary layer. 



5.1.2. Algorithm [2} We now consider Algorithm [2] (analyzed in Section 4.3), 
where the reconstruction at the end of each parareal iteration is performed using the 
matching operator V. 

The maximal number of parareal iterations is set at K — 6. Figure 5.2 shows 
the macroscopic and microscopic errors as a function of e, for the chosen values of 
k. The numerical results are in agreement with Theorem |13[ At each odd parareal 
iteration, the order of convergence (in terms of e) of the macroscopic error increases 
by 1, whereas the microscopic error decreases, but remains of the same order in e. 
At each even iteration, the converse holds: the order of convergence (in terms of e) 
increases by 1 for the microscopic error, whereas the macroscopic error decreases but 
its order remains alike. Note also that, for the smallest considered values of e, the 
algorithm reaches machine precision in 5 to 6 iterations. 

As with Algorithm [T] when e is too large, the numerical results do not agree 
with the theoretical results, because the chosen time-step At does not satisfy the 
assumption At > tf L . 



At this point, we have numerically verified our theoretical results, and know that 
the macroscopic error is bounded from above by, and actually roughly of the order of, 



sup|££| 

n 



At (id 
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-ffc/21 



(5.2) 




where Ck,At a priori depends on k and At, but is independent of e (and likewise for 
the microscopic error). 

On Figure |5.3| we plot the macroscopic and microscopic errors as a function of 
the iteration number k, < k < K = 30, for various values of e. We observe an 
exponential convergence to the exact solution as a function of k, with a convergence 



rate that increases when e decreases. We deduce from (5.2) how the constant Ck At 



depends on k: there exists C\ t and Cj± t , independent of k and e, such that 



supk: 



'Pit) 



-r/s/21 / 

\AtJ 



(5.3) 



Note that, for e = 10 1 (which is quite a large value compared to At — 10 1 ), 
convergence is very slow, which is in agreement with the previous observations. 
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Fig. 5.3. A Igorithm \^}or the system l |5.1| l, with exact fine-scale and coarse propagators and 
parareal time step At = 10" 1 / errors as a function of k for different values of e (left: macroscopic 
error; right: microscopic error). 



To understand how C\ t and C\ t depend on At, we perform another experiment, 
in which we fix e = 10~ 5 and vary At. We then plot the error as a function of At^ 1 for 



different values of k (see Figure 5.4). These results show that the macroscopic error 
varies proportionally to (At -1 ) + Combined with (5.3), we therefore deduce 
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macro 



that, on this numerical test-case, the macroscopic error satisfies 

SU P \Ek I ~ Cmacro ( Cm 
n 

and likewise for the microscopic error: 

SU P l e fcl ~ ^micro (^micro ^) 



l+Lfc/21 



(5.4) 



(5.5) 



We in particular see that, if e/ At is sufficiently small, then the parareal trajectory 
converges to the exact trajectory when k goes to oo. 
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Fig. 5.4. Algorithm^for the system \5.1\ , with exact fine-scale and coarse propagators and 
e = 10~ 5 : errors as a function of At -1 for different values of k (left: macroscopic error; right: 
microscopic error). 



5.1.3. Algorithm 3. To complete these numerical tests, we consider Algorithm 



3 originally proposed in [51 [2H] (which we analyzed in Section 4.4), and repeat the 
previous experiment. The results, shown in Figure |5.5| are in agreement with the 
theoretical results. 
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Fig. 5.5. Algorithm 3 for the system with exact fine-scale and coarse propagators: errors 

as a function of e for different values of k (left: macroscopic error; right: microscopic error). 
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5.2. Results with exact microscopic and approximate macroscopic in- 
tegration. In this section, we repeat the above experiments, but now using a forward 
Euler time discretization for the coarse propagator, using the time step At (hence, 
to propagate the system over the time range At, the scheme Ca* consists in doing a 
single step of the forward Euler algorithm). The fine-scale propagator is again the 
exact one. 

5.2.1. Algorithm [T| We first consider Algorithm [I] (for which the reconstruc- 
tion is done using the lifting operator C), and set the maximal number of parareal 
iterations to K = 3. Figure 5.6 shows the errors as a function of e for the chosen 
values of k. 
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FlG. 5.6. Algorithm\l\for the system j5.1) , with exact fine-scale propagator and approximate 
coarse propagator: errors as a function of e for different values of k (left: macroscopic error; right: 
microscopic error). 



When comparing Figure 5.6 with Figure |5T| (in which the macroscopic dynamics 



is exactly integrated), we notice that the behavior of the algorithm is similar for large 
values of e. For small values of e, the errors approach an asymptotic value as e goes 



to zero (rather than converging to as in Figure 5.1 1, with an asymptotic value that 



depends on the number of iterations k. The larger k is, the smaller this asymptotic 



value is, and the wider the range of e for which results of Figures 5.6 and 5.1 (with 
approximate or exact integration at the macroscopic scale) agree. 



This observation is confirmed in Figure 5.7 where we show the errors as a function 
of the iteration number k, 1 < k < K = 8, for various values of e. We see that the 
errors first converge exponentially to as k increases, and then reach a plateau. The 
residual macroscopic (resp. microscopic) error is of the order of 0(e 2 ) (resp. 0(e)). 

We explain this behavior as follows. The parareal iterations iteratively correct 
the approximation made using the coarse propagator. When the coarse propagator is 
a forward Euler discretization of the approximate macroscopic equation, there are two 



sources of error: a modeling error (due to the fact that the macroscopic equation (2.4 1 
is only an approximation of the slow part of the reference model (2.1 )), and a time 
discretization error. For large values of e, the modeling error is dominant, and the error 
behaves as if the coarse propagator were exact. For small e, the time discretization 
is dominant, and the error becomes therefore independent of e. Due to the parareal 
iterations, the time discretization error is iteratively removed. However, due to the 
fact that the reconstruction is performed using the lifting operator C, the modeling 
error never vanishes when e > 0. Hence, when k goes to infinity, Algorithm [T] using an 
approximate coarse propagator converges to the solution given by a parareal algorithm 
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FlG. 5.7. Algorithm^for the system with exact fine-scale propagator and approximate 

coarse propagator: errors as a function of k for different values of e (left: macroscopic error; right: 
microscopic error). 



with no time-step discretisation error (this latter has been removed by the iterations 
in k), but with some modeling error. The solution hence converges to that given by 
Algorithm [T] with exact coarse propagation. 

5.2.2. Algorithm [2} We now consider Algorithm [2] (for which the reconstruc- 
tion is performed using a matching operator V) 1 and set the maximal number of 
parareal iterations at K — 6. Figure |5.8| shows the errors as a function of e for the 
chosen values of k. 
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Fig. 5.8. Algorithm^ for the system \5.1\ , with exact fine-scale propagator and approximate 
coarse propagator: errors as a function of e for different values of k (left: macroscopic error; right: 
microscopic error). 



We compare Figure 5.8 to the corresponding Figure 5.2 (where the macroscopic 
equation is exactly integrated). We again notice that, for small e, the algorithm be- 
haves differently: in particular, the convergence when e goes to zero slows down when 
the macroscopic equation is only approximately integrated. However, the behavior 
with respect to k is left unchanged. We show on Figure 5.9 the evolution of the errors 
as a function of the parareal iteration number k, < k < K = 30, for a number of 
fixed values of e. As in Figure |5.3| the computed trajectory again converges to the 
exact microscopic solution up to machine precision, exponentially with respect to k, 
despite the presence of time discretization errors at the macroscopic level. Moreover, 
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comparing these results with those obtained when using an exact coarse propagator 
(see Figure 5.3 1, we see that only a few extra parareal iterations are needed. 




10 20 30 10 20 30 

k k 

Fig. 5.9. Algorithm^ for the system with exact fine-scale propagator and approximate 

coarse propagator: errors as a function of k for different values of e (left: macroscopic error; right: 
microscopic error). 



6. Nonlinear examples. We finally illustrate the performance of our Algo- 
rithm [2] on two nonlinear examples. Such cases are not covered by the theoretical 
analysis of Section |4j where we considered a linear problem. 

The first nonlinear example we consider is a straightforward extension of prob- 



lem (2.1 1, and reads 



= -Aa 



V, V= ~{x - y), 



which is of the form (2.5 1. The corresponding macroscopic model is 



X 



-XX - x 1 



(6.1) 



(6.2) 



We use Algorith m [2| t o integrate this system. The fine propagator .Fa* is a forward 
Euler scheme for with the time step 5t = 10~ 5 . The coarse propagator Ca* is a 
forward Euler scheme for (6.2 1 with the time step At = 10 _1 (which is equal to the 
parareal time step). The lifting operator reads C(X) — (X,X 2 ), and the matching 
operator is again given by (3.9 1. The remaining parameters are chosen identical to 
those in the previous experiments. On Figure |6.1[ we show the error as a function 
of e for different values of k. Comparing that figure with Figure 5.2 we see that 



Algorithm [2] performs equally well on this nonlinear case as on the linear problem 
considered in Section 

The second nonlinear example we consider is the so-called Brusselator problem, 
which was already considered in e.g. |19j . It reads 



y 



yx\ 



x\x 2 . 



(6.3) 



-(B -y) 

e 



yx x . 



It models the evolution of the concentration of three chemical species. The concen- 
tration of y is reduced via reaction with x, but restored to its base level Bq with a 
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function of e for different values of k (left: macroscopic error; right: microscopic error). 



characteristic time of the order of e. We choose A = 1 and Bq = 3. The fine propaga- 
tor T^t is a forward Euler discretization of (6.3) with the time step St — 10~ 5 . The 



coarse propagator Ca* is a forward Euler discretization of the macroscopic model 



±x = A - (B + l)asi 



x\x 2 , 



X 2 = Bq Xi 



x\x 2 , 



with the time step At = 1CP 1 (equal to the parareal time step). This system thus 
has two slow variables and one fast one: u = (x,y), with x — (x\,X2)- Note that 
this case does not enter our theoretical framework for two reasons: (i) the example 
is nonlinear; and (ii) the equation for y in the microscopic model is not purely a fast 
equation (the second term in the right-hand side of the equation for y in (6.3) is not 
scaled by e _1 ). 



We show on Figure 6.2 the results obtained. Algorithm[2]again performs very well. 
Actually, on this problem, the convergence behavior of Algorithm [2] closely resembles 
that of Algorithm 3: at parareal iteration k, the order of convergence (in terms of e) 
seems to be equal to k, both for the macroscopic and the microscopic errors. 



— 10 




Fig. 6.2. A Igorithm p| for the system l |6.3| , with parareal time step At = 10 1 : errors as a 
function of e for different values of k (left: macroscopic error; right: microscopic error). 



7. Discussion and conclusions. We have introduced and analyzed two micro- 
macro parareal algorithms for the time-parallel integration of singularly perturbed 
ordinary differential equations, and provided a numerical analysis for the case where 
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the problem is linear and the coarse and fine-scale propagators integrate the macro- 
scopic, resp. microscopic models exactly. The analysis shows that, when an appropri- 
ate matching operator is used to update the microscopic state after correction of the 
macroscopic state (which corresponds to using Algorithm [2]) , the rate of convergence 
(in terms of the modelling error e, which quantifies the time scale separation between 
the microscopic and the macroscopic evolutions) improves at each parareal iteration, 
and is roughly equal to e k l 2 . We have illustrated this theoretical result with numerical 
experiments, and also numerically investigated the effect of using a numerical scheme 
to integrate the macroscopic model (thereby introducing some discretization error). 
The results show that the proposed micro-macro parareal algorithm, Algorithm [2j 
is robust with respect to time discretization errors at the macroscopic level. It can 
hence be viewed as an interesting way of using an approximate macroscopic model to 
speed up simulations of high-dimensional multiscale systems of singularly perturbed 
ordinary differential equations, provided that special care is taken when transferring 
information between the microscopic and macroscopic levels of description. 

Several questions remain open. First, while the analysis reveals that it is impor- 
tant to choose the parareal time step At sufficiently large (to average out the initial 
time boundary layers in the full microscopic dynamics) , the dependence on At of the 
numerical error and the convergence rate have not been analyzed. In particular, one 
may expect an optimal time step At to exist that leads to a required accuracy with 
a minimal cost. Assume again (as for the original parareal algorithm, see the intro- 
duction) that the cost of a single evaluation of the fine-scale propagator is much 
larger than the cost of propagating the macroscopic system using C& t over the com- 
plete time range [0, T] (This assumption is all the more justified as the macroscopic 
system is a low-dimensional problem compared to the microscopic problem). Then 
the cost of the parareal algorithm, after K iterations, is proportional to KAt/e (We 
have assumed that the cost of J- At is proportional to At/e, since we need to use a time 
step of the order of e over a time range of length At). This cost is to be compared 
with the cost of the full microscopic sequential integration, which is proportional to 



NAt/e. The computational speed-up is thus N/K. We saw on Figure 5.2 that, for 
reasonably small values of e, results obtained at the iteration K — 6 were satisfactory. 
For the test-case considered in Section |5.f .2[ the computational speed-up is thus 

N Tl At 

k = J iT = 16 - 6 ' 

Second, we expect Algorithm [2] to extend to more general dissipative systems. As 



pointed out above, we considered here the simple linear problem (2.1 1 to focus on the 
issues stemming from using two different levels of description of the same system. We 
have already checked in Section [6] that Algorithm [2] behaves equally well on nonlinear 
systems of singularly perturbed ODEs. We currently investigate the extension of the 
algorithm to a setting, motivated by molecular simulations, where the reference model 
is a high-dimensional stochastic differential equation (modeling the evolution of all 
the degrees of freedom of the atomistic system) and the macroscopic model is the 
effective dynamics of the slow component of the microscopic model, derived under 
time scale separation assumptions following 24, 25] , 
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Appendix A. Proofs of Lemma [2] and Corollary [3} 

Before proving Lemma[2]and Corollary[3j we start with a preliminary result. Here 
and in all what follows, || • || denotes the Euclidean norm when applied to vectors, and 
the associated operator norm when applied to matrices. 
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Lemma 15. Let M be a matrix in R dxd such that the real part of the spectrum 
of M is positive. Then, there exist C > and /i > such that, for all time t > 0. 

|| exp(-Af£)|| < Cexp(-/xt). (A.l) 

One can choose /i = inf„ g(7 (jw) R- e (^)/2, where a(M) denotes the spectrum of M. We 
also have 

\\M^\\ <-■ (A.2) 
(j, 



Proof. We introduce the Jordan form of the matrix M. Let us assume for simplic- 
ity of notation that M has only two Jordan blocks associated to two complex eigenval- 
ues Ai 7^ A2 with < Re(Ai) < Re(A2). The generalization to any number of Jordan 
blocks is straightforward. Let us denote fi = vai^^ a (M) Re(f)/2 = Re(Ai)/2 > 0. 

Since M has only two Jordan blocks, there exists an invertible matrix Q € M. dxd 
such that 



M = Q- 



N dl (\x) 
N d2 (X 2 ) 



Q. 



where d\ + d 2 = d, and, for any m € N* and any A G 



N m (X) 



A 1 
A 



A 1 
A 



We compute, for any ieR, 

exp(Mt) = Q- 1 



exp(Aii)P dl (t) 







exp(A 2 *)P d2 (t) 



Q, 



(A.3) 



where, for any m € N* and any t G K, 



1 f t 2 /2 t 3 /6 
1 i t 2 /2 



t m - l /{{m- 1)!) 
t m - 2 /((m - 2)!) 



Since P m (t) is a matrix with entries which are polynomial functions of f, there exists 
a constant C that only depends on the matrix M such that 

Vt > 0, ||P dl (-t)|| + ||Pd,(-t)|| < Cexp(jjt). 



We then infer from (A.3) that there exists a constant that only depends on M such 
that 



Vt > 0, ||exp(-Mt)|| < Cexp(-fd). 
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This yields (A.l). Then, (A. 2 1 is obtained using the fact that 



Im- 1 !! = 



poo poo 

/ exp(-Mt)dt < / ||exp 
Jo Jo 



C°° C 
>(— Aft)|| dt<C exp(-fd)dt = -. 

Jo A* 



This concludes the proof of Lemma [15] □ 

We are now in position to prove Lemma [2] and Corollary [3] 
Proof of Lemma [2} We start by writing 

i = y — A~ 1 qx 
— - (qx — Ay) — A~ x q {ax + p T y) 



A 



- + (A- 1 q)p 1 



z-X(A- 1 q). 



(A.4) 



where A is defined by (2.4). Introducing 

M e := A + e (A- X q) p T e IR(<*-i)x (<*-!), y ■- \ (A^q) G 



we recast (A.4) as 



M e 



-z — Vx. 



(A.5) 



From the definition of M e , and in view of Assumption ( |2.3[ ), it is clear that there 
exists a critical value eo(A, q,p) such that for all e < eq(A, q,p), the matrix M e has a 
spectrum with a real part bounded from below by A_/2 > 0, where A_ is independent 
of e. Up to a modification of 6q(A, q,p, a), the same property holds true for the matrix 
M e + eAId that will appear below (where A is defined by ( |2.4[ )). In the sequel of the 
proof, we will systematically work with e < eo(^4, q,p, a). 



By explicit integration of (A.5), we have 



z(t) - exp(-MH/e)z 



cxp [-M € (t - s) /e] V x(s) ds. 



(A.6) 



From (2.1 ), we have x — \x + p z. Using equation (A.6), we thus obtain 



x(t) — xo exp(Xt) = p / exp (X(t — s)) z(s)ds 
Jo 

— p T exp (X(t — s)) exp (— M e s/e) zods 
Jo 

-p T ( exp(A(t-s)) / exp(-M e (s-r)/e)Vx(r)drds. (A.7) 



(i 
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To bound the first term of (A. 7), we write, using Lemma 15 



= exp(At) 

< 

< e 

< £ 



[ exp(A(i- s))exp(-M e s/e)ds 
Jo 

f exp (— (M e /e + Aid) s) ds 
Jo 

(M e /e + Aid) -1 ||exp(At)Id-exp(- [M e /e)t)\\ 

(M e + eAId) _1 (||exp(Ai)Id|| + ||exp (- (M e /e) t)\\ ) 
C{A ' q ^ a) ( exp(AT) + C(A, q, p, a) exp (-A_t/(4e)) 
< C{A,q,p,a,T)e, 



(A.8) 



when e < eo(A, q,p, a). Turning to the second term of (A. 7 1, we use Fubini's theorem, 
and write 



exp(A(t-s)) / exp(-M e (s- r) / e) V x(r) drds 



exp(At) / exp(M e r/e) 



/ exp(-(M £ /e + \ld)s)di 

J r 



V x(r) dr 



exp(Ai) / exp(M e r/e) exp (-(M e /e + AId)r) - exp (-(M e /e + Xld)t) 
Jo L 

x {M f - / e + Xl&y 1 V x{r) dr 
[exp (A(i - r)) Id - exp (~M e (t - r)/e)] (M e /e + Aid) -1 Va;(r) dr. 



Therefore, for e < eo(A, q,p, a), using Lemma 15 we obtain 



exp(A(t — s)) / exp (— M e (s — r)/e) Vx{r)drds 



< 



< e 



(MVe + AId)" 1 II^H sup \x(r)\ \\exp(\(t - r))ld- exp(-M £ (t - r)/e)\\dr 

0<r<t Jo 



(Ar + eAId) 1 C(A,q,p,a,T) sup \x(r) 

0<r<t 



<C(A,q,p,a,T)e sup \x(r)\. 

0<r<t 



(A.9) 



Combining equations (A.7|, (A.8I and (A.9), we get 



\x(t) - x exp(At)| < C(A,q,p,a,T)\\p\\ \\z \\ e + C(A,q,p,a,T) \\p\\ e sup \x(r) 

0<r<t 

<C{A,q,p,a,T)e(\\z Q \\+ sup \x(r)\ J , 

\ 0<r<t I 
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and hence, 



sup \x(t) — xq exp(At)| 

0<t<T 

< C(A,q,p,a,T) e[\\zo\\+ sup \x(r)\ 

\ 0<r<T 



< C(A, q,p,a,T) e [ \\z \\ + \x \ + sup \x{r) - x exp(Ar)| 

V 0<r<T 

We deduce that, for c < co(A, q,p, a, T), 

sup \x(t) -xoexp(Xt) \< - — — — (\\z + Fo ) 

o<t<T 1 - C(A,q,p,a,T)e 

<C{A,q,p,a,T) c (||2 || + |a;o|). 



This proves {2.8 1. 

We now turn to proving the bound ( |2.9[ ) on z(t). Introducing 

B := (A- 1 q)p T eK («i-i)x(<i-i) ) 



we now recast (A.4| as 



.4 



D 



z - Vx = z - [Bz + Vx] 

e 



By explicit integration, we have 



z(t) - exp(~At/c)z = cxp [-A(t - s)/e] (Bz(s) + Vx(s)) ds. 



Using (A.l) and Assumption (2.3|, we obtain 



\\z(t) - exp(-At/e)z \\ < C{A) / exp 



s - t 



2c 



< C{A,q,p,a) / exp 



\\Bz(s) + Vx(s)\\ds 
s - t 



2c 



(A.10) 



(A.ll) 



<C(A,q,p,a) sup (||z(s)|| + |a;(a)|) 
se[o,t] 



(\\z(s)\\ + \x(s)\)ds 
2c 



< cC(A, q,p, a) I sup \\z(s) — exp(— As/e),zo|| 
\se[o,t] 



+ sup || exp(— As/c)zq\\ + sup \x(s) 
se[a,t] se[a,t] 

Taking the supremum over t £E [0, T], we obtain, for e < cq(A, q,p, a), 

eC(A,q,p,a) 



sup \\z(t) -exp(-At/e)z \\ < 
te[o,T] 1 - tC{A,q,p,a) 



< 



eC{A,q,p, a) 
1 - cC(A,q,p,a) 
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sup || exp(— ^4s/e)zo|| + sup \x( 

\s£[0,T] sS[0,T] 

sup \x(s)\ ) . 



We then deduce from (A.10) that, for e < e (A,q,p,a,T), 

sup \\z{t) - exp(-At/e)z < — 

te[o,T] 1 - eC(A,q,p,a 



N) 



<eC(i4,ff,p ) o,T)(||z || + |a:o|). 



(A.12) 



This proves (2.9 1. 



We finally turn to proving (2.11). Using (A.l), we see that 
||exp(-A*/e)|| < C(A)exp(-A_t/(2 e )), 



2c 



thus, for times t > tf = — ln(l/e), we have || exp(— .Af/e)|| < C(A)e. We then 
deduce from (A.12) the bound (2.11). This concludes the proof of Lemma|2] 

Proof of Corollary [3} The first assertion follows directly from ( |2.8[ ) and the 

fact that ||zo|| < \\yo\\ + C|a;o|. The second assertion follows from (2.11) and (2.13). 
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